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We assess the potential of two-terminal graphene-hBN-graphene resonant tunneling diodes as 
high-frequency oscillators, using self-consistent quantum transport and electrostatic simulations to 
determine the time-dependent response of the diodes in a resonant circuit. We quantify how the 
frequency and power of the current oscillations depend on the diode and circuit parameters including 
the doping of the graphene electrodes, device geometry, alignment of the graphene lattices, and the 
circuit impedances. Our results indicate that current oscillations with frequencies of up to several 
hundred GHz should be achievable. 


Resonant tunneling diodes (RTDs) operating at fre¬ 
quencies and output powers of up to 1.4 THz and 10 pW 
have been recently demonstrated [IH3] . A new addition to 
the family of devices that exhibit resonant tunneling and 
negative differential conductance (NDC) is the graphene- 
based tunnel transistor [iHT5]. In this device, a thin sheet 
of hexagonal boron nitride (hBN) acts as a potential bar¬ 
rier separating two monolayer graphene electrodes. The 
NDG arises from the constraints imposed by energy and 
momentum conservation of Dirac fermions, which tunnel 
through the hBN barrier when a bias voltage is applied 
between the graphene electrodes. Peak-to-valley current 
ratios (PVRs) of 2:1 have been seen at room tempera¬ 
ture, with peak current densities of 0.28 pA/pm^ [HI [7]. 
When these devices are placed in a resonant circuit and 
biased in the NDC region, MHz oscillations occur [7]. 

Here, we use a theoretical analysis to investigate how 
the device and circuit parameters can be tuned to in¬ 
crease the operating frequency of graphene resonant 
tunneling diodes (GRTDs). Our model device, shown 
schematically in Fig. 1(a), comprises two graphene lay¬ 
ers separated by a hBN tunnel barrier of thickness, d. 
The bottom (B) and top (T) graphene electrodes are ar¬ 
ranged in an overlapping cross formation, resulting in an 
active tunneling region of area, A = 1 pm^. We con¬ 
sider the general case when the two graphene crystalline 
lattices are slightly misorientated by a twist angle, 6, 
see Fig. 1(a). The resonant tunnel current is particu¬ 
larly sensitive to this angle [7j. A bias voltage, 14, ap¬ 
plied between the top and bottom graphene layers [Fig. 
1(b)] induces a charge density, pb,t, in each layer and 
causes a tunnel current, /&, to flow through the hBN 
barrier. The graphene layers, with in-plane sheet resis¬ 
tance, R, carry current, 7, (black arrows) from two pairs 
of Ohmic contacts [orange in Fig. 1(a)] to the central 
active (tunneling) region of the device, i.e. currents, //2, 
flow to/from each contact. The electrostatics of the diode 
[1] are governed by the equation et4 = Pb — Pt + <Pb, 
where (jjb = eF^d is the electrostatic potential energy dif¬ 
ference across the barrier, with 74 being the electric field 
in the barrier, e is the magnitude of the electronic charge, 
and pb,t are the two Fermi levels [see Fig. 



FIG. 1. (a) Schematic diagram of the GRTD comprising bot¬ 
tom (red) and top (blue) graphene lattices, misaligned by an 
angle 6 and separated by a hBN tunnel barrier (dark green). 
The current, 7, passes through the tunnel barrier between the 
graphene electrode layers to/from Ohmic contacts (orange). 
The diode is mounted on a hBN layer (light green) and an 
insulating substrate (purple), (b) Schematic diagram of the 
resonant circuit incorporating the GRTD (in box) showing 
the voltage applied, V, circuit inductance, L, and resistance, 
R. The band diagram is shown (box), with electrostatic pa¬ 
rameters defined in the text. 

A device with NDG can generate self-sustained current 
oscillations when placed in an RLG circuit m- To inves¬ 
tigate the frequency response of the GRTD, we solve the 
time-dependent current continuity and Poisson equations 
self-consistently, using the Bardeen transfer Hamiltonian 
method to calculate the tunnel current, 

E \MWB{EB)-fT{ET)]SiEB-ET-M, 

ks -Mt 

( 1 ) 

as a function of time, t, and 14. The summation is 
over all initial and final states, with wavevectors, \i.B,Ti 
measured relative to the position of the nearest Dirac 
point in the bottom layer, = (±47r/3ao, 0), where 
± distinguishes the two non-equivalent Dirac points in 
the Brillouin zone and oq = 2.46 A is the graphene lat¬ 
tice constant. The Fermi function in each electrode is 
IbHEb^t) = [1 + where Eb,t = 

SBprhvpkBpr is the electron energy and sb,t = ±1 labels 
electrons in the conduction (-I-) and valence (—) bands, 
at temperature T = 300 K. Tunneling between equiv¬ 
alent valleys gives the same contribution to the tunnel 
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current, so we consider transitions between K+ points 
only. In Eq. Q the matrix element, M, is 

M = S7(6»)5((pB, ^T)Vs{ci - AK), (2) 

where S is a normalisation constant, 7(0) is the spatial 
overlap integral of the cell-periodic part of the wavefunc- 
tion, g{^PB, ‘Pt) describes electron chirality. Vs is the elas¬ 
tic scattering potential, and q = — ky (see below). 

In recently-studied GRTDs [7], the crystalline lattices 
of the two graphene layers are misorientated by only a 
small twist angle 0 « 1°. Nevertheless, this gives rise 
to a significant misalignment of the Dirac cones of the 
two layers, AK = {R{0) — 1)K'^, where R{d) is the 2D 
rotation matrix. When 9 < 2°, |q| « |AK| = AK and 
electrons tunnel with conservation of in-plane momen¬ 
tum. However, tunneling electrons can scatter elastically 
from impurities and defects, broadening the features in 
the current-voltage curves. Therefore, we use a scatter¬ 
ing potential Vs{q) = Vo/{q‘^ -I- qc^), with characteristic 
lengthscale 1/qc = 15 nm, which gives the best fit in the 
region of the resonant peak and NDC [5]. 

The misorientation of the layers also causes the cell- 
periodic part of the Bloch wavefunctions in the two layers 
to become misaligned, thus reducing their spatial over¬ 
lap integral, 7(0), and, consequently, the tunnel current 
amplitude. We estimate 7(0) by calculating the overlap 
integral of the normalised cell-periodic parts of the Bloch 
states, u{r), at the Dirac point in the two electrodes over 
an area. Sc, that greatly exceeds the length scale, ^ l/qc, 
of the impurity potential, i.e. we take 

7(0) =/ dScu*(R{9)r)u{r). (3) 

Jsc 

The chiral wavefunctions give rise to the term 

9{Pb,Pt) = l + SBe^^^+STe-^^^+SBSTe^^‘^^-^^\ ( 4 ) 

where (p = {ky/kx) is the orientation of the 

wavevector. Finally, in Eq. ([^, S = where ^ 

depends on the wavefunction amplitude in the barrier, 
K = y/2mAi, jh, is the decay constant of the wavefunction 
in the barrier, of height Af, = 1.5 eV, and m = O.Srrie is 
the effective electronic mass in the barrier [Ij. 

Fig. i shows the equilibrium (static) current-voltage 
hiVk) curve (blue), where Vf, ^ V and Ih = I, calcu¬ 
lated for an undoped device with 6 = 0.9° and a barrier 
width d = 1.3 nm (4 layers of hBN). These parameters 
are similar to those of a device which has been recently 
fabricated and measured [7]. The calculated It (14) curve 
reproduces the line-shape, position of the resonant peak 
and current amplitude of the measured device charac¬ 
teristics ill]. The peak in current occurs when many 
electrons can tunnel with conservation of momentum, i.e. 
q— AK ^ 0, corresponding to a resonant increase in the 
matrix element M, i.e. when = hvpAK, for 9 close 
to 1°. Our simulations and the measurements ill! indi¬ 
cate that temperature has negligible effect on the current- 
voltage curve when I4 > kT/e ^ 30 mV. Phonon-assisted 



FIG. 2. Equilibrium and non-equilibrium current-voltage 
curves calculated for 6 = 0.9°, L = 140 nH and R = 50 Q. 
Blue curve: equilibrium current-voltage characteristic IbiVb). 
Note, in equilibrium, Vb ~ V and h = I- Green curve: time- 
averaged current (/(t))t vs V. Red curve: peak-to-peak volt¬ 
age amplitude (right scale) of the stable current oscillations. 
Inset: I{t) plot showing stable oscillations with / = 4.2 GHz. 


tunneling has a relatively weak effect on the device char¬ 
acteristics when the graphene crystal lattices are almost 
aligned [7j. 

We now consider the non-equilibrium charge dynamics 
when the device is in a series circuit with inductance, L, 
and resistance, R, see Fig. 0b) ; the diode has its own 
in-built capacitance, C. The primary contribution to R 
arises from the graphene electrodes |3] and depends on 
the charge densities, pb,t- This dependence does not 
have a significant effect on the high-frequency (HF) re¬ 
sponse: for most of the oscillation period, changes in pb,t 
do not greatly affect R. Therefore, to reasonable accu¬ 
racy and for simplicity, we take R to be independent of 
time. However, the value of R can be changed by altering 
the device geometry, for example, by reducing the length 
of the electrodes, and we consider this effect on the per¬ 
formance of the GRTD. We also consider how L affects 
the oscillation frequency, which could be controlled by 
careful design of the microwave circuit, for example by 
using a resonant cavity or integrated patch antennas [5] . 

We determine the current, I{t), in the contacts and ex¬ 
ternal circuit by solving [T7j self-consistently the current- 
continuity equations: dpB,T/dt = ±(/b — I)/A, where 
the -I- (—) sign is for the bottom (top) graphene layers, 
see Fig. 1(b), pb,t are related by Poisson’s equation: 
ePb = Pb - Pbd = -{pt - Ptd), in which e = eoCr- and 
er = 3.9 H Hg is the permittivity of the barrier, and 
Pbd (ptd) are the doping densities in each layer. The 
voltages across the inductor and resistor, V4 and Vr, 
and the currents through them, Ib and /r, are given by 
dl/dt = Vl/L, Vr = IR, and V = Vr + V, + Vr. 

Following initial transient behavior, I{t) either decays 
to a constant value or oscillates with a frequency, /, and 
time-averaged current, Fig. inset, shows a 

typical I{t) curve, for V = 0.48 V, exhibiting current os¬ 
cillations with / = 4.2 GHz. In Fig. we show (/(t))t 
versus V (green curve) and IbiVb) characteristics (blue 
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FIG. 3. (a) fmax{R) calculated when Nl = 4. Inset: Log-log 
plot, (b) fmax vs R when Nl = 2 (red), 3 (green) and 4 
(blue). Inset: fmax vs Nl calculated when R — 50 ^1. Curves 
are shown solid over the range of R presently obtainable in 
GRTDs and dashed for R values that could be achieved by 
future device designs. All curves are for undoped devices. 

curve) for an undoped device, with 6 = 0.9°, placed in 
a resonant circuit with i? = 50 11 and L = 140 nH. The 
plot reveals that when V is tuned in the NDC region 
(0.55 V < R < 0.8 V), AI4 = (red curve) 

becomes non-zero indicating that self-sustained current 
oscillations are induced. Here, jg niaxi- 

mum/minimum voltage dropped across the inductor dur¬ 
ing a current oscillation period. Also, the {I{t))t versus 

V curve (green) diverges from the equilibrium current, 
h{Vb) (blue curve) in the NDC region. This is due to 
asymmetric rectification of I(t) in the strongly nonlin¬ 
ear NDC region of Ib{Vb). When the device is biased in 
regions of positive differential conductance (PDC), i.e. 

V < 0.55 V or H > 0.8 V, the current oscillations are 
suppressed and (/(t))t converges to Ib{Vb)- 

This behavior is similar to that recently measured in 
a GRTD, where current oscillations with / ~ 2 MHz 
were reported[7j. That device had high circuit capaci¬ 
tance due to the large-area contact pads and coupling to 
the doped Si substrate (gate). This effect can be mod¬ 
elled by placing a capacitor in parallel with the GRTD. 
Including this large capacitance (65 pF) limits the max¬ 
imum frequency, as observed [7]. Here, we consider the 
case when parasitic circuit capacitances are minimised, 
using the geometry exemplified in Fig. 1(a), so that the 
only significant contribution to the total capacitance is 
from the graphene electrodes, as described by the charge- 
continuity equation. This enables us to investigate the 
potential of GRTDs optimised for HF applications. 

A small signal analysis[T2] provides insight into how L, 
i?, and the form of Ib{Vb) affect the circuit response and 
gives an approximate oscillation frequency: 

r = h\l(1 - r/Rn) - q-n (1 - qir/Rn? /4, (5) 

where Rat is the maximum negative differential resistance 
of the equilibrium I{V) curve, the circuit factor Qjq = 
Rn\/CJL, and /o = 1I2tt^/LC. For our device, Rn is 
large and therefore ~ /o- For a given C (that depends 


on A and d), the oscillation frequency can be increased 
by reducing L. The decay parameter of the small signal 
analysis reveals that the circuit will oscillate only if 

{RnIR-QI)>Q. ( 6 ) 

Consequently, R, and the shape of the static Ib{Vb) curve 
are also important for optimising the HF performance. 

We now consider the fully self-consistent simulation of 
the charge dynamics obtained using Eqs. (1-4). Fig. 
[^a) shows the fmax{R) curve calculated for the diode 
parameters, which compare well to recent measurements 
[7], used to produce the /&(t4) curves in Fig. We 
determine fmax{R) by finding the smallest L value for 
self-sustained current oscillations. The solid part of the 
curve in Fig. 3(a) shows fmax over the range of R values 
that can be achieved by only small modifications to the 
design of existing devices, for example by reducing the 
length of the graphene between the tunnel area and the 
Ohmic contacts, or doping the electrodes. The dashed 
part of the curve is calculated for R values that may be 
possible in future configurations. The curve reveals that 
for a readily-attainable i? = 50 H, fmax = 1-8 GHz. 

Fig. 3(a), inset, reveals the power law fmax cx 
which can be derived by setting Eq. (|^ equal to zero 
and rearranging to find the smallest L value for a given 
R, Rat, and C [16]. For this case 

= (27rC/^)-i (X R-°■^ (7) 

which compares well with the full signal analysis. 

To increase fmax, in addition to varying the external 
circuit parameters, we can also modify the Ib{Vb) curve. 
Reducing the number of hBN layers, Nl, in the tunnel 
barrier significantly increases the tunnel current (~ 20x 
for each layer removed 119j) thus reducing Rat and in¬ 
creasing fmax, see Eq. ([^Fig. [^b) shows fmax{R) cal¬ 
culated for a device with Nl = 4 (blue curve), 3 (green 
curve) and 2 (red curve). Reducing the barrier width 
produces a large gain in fmax for all R. For example, 
fmax for a device with Nl = 2 is at least an order of 
magnitude higher than when Nl = 4 (e.g. for R = 50 fl, 
fmax= 26 GHz when Nl = 2, compared to fmax = 1-8 
GHz when Nl = 4). 

The Ib{Vb) characteristics can also be modified by dop¬ 
ing the graphene chemically [SniET] or, equivalently, by 
applying a gate voltage, Vg, to shift the position of the 
current peak and, thereby, strongly influence Rat and 
the peak to valley ratio laiT]. However, a gate electrode 
would capacitively couple to the graphene layers and the 
contact pads, so we do not consider it here. In Fig. |^a), 
we show Ib{Vb) curves calculated when Nl = 2 for an 
undoped (red curve) and an asymmetrically-doped de¬ 
vice with pBoje. = 10^^ cm“^ and prole. = 0 (green 
curve). When pbd > 0, the resonant peak occurs at 
higher 14 than when pbd =0, and the magnitude of the 
current peak is also higher, raising the PVR from 1.5 to 
3.5. 

The shoulder of the green curve in Fig. 4(a), indi¬ 
cated by an arrow when pBo/e = 10^^ cm“^, arises from 
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FIG. 4. (a) Ib{Vb) characteristics calculated for a doped (green 
curve, Pbd/& = 10^^ cm“^) and undoped (red curve) device, 
with Nl = 2. The arrow shows the shoulder that arises due 
to the quantum capacitance effect, (b) fmax vs R curves for 
the devices in (a). Inset: fmax vs [pboI^) calculated when 
i? = 50 fl, with ptdIc = 0. Curves in (b) are shown solid over 
the range of R presently obtainable in GRTDs and dashed for 
R values that could be achieved by future device designs. 


the low density of states around the Dirac point. This 
gives rise to an additional quantum capacitance [SI 
Cq, whose effect is most prominent when the chemical 
potential in one layer aligns with the Dirac point in the 
other layer. The total device capacitance is given by 
C~^ = Cq^ +Cq^, where Cq = t^trAjd is the geometric 
capacitance. When p,B,T passes through the Dirac point, 
Cq -p- 0 and, hence, C —)■ 0, suggesting that the RC time 
constant of the device could be reduced. In practice, Cq 
is small for only a small fraction of the oscillation period 
and so its effect on the fundamental frequency of I{t) is 
negligible. 

Fig. I^b) shows the fmax{R) curves calculated for the 
undoped (red) and doped (green) devices and reveals that 
the doped device is significantly faster for all R. Fig. 
I^b) inset shows that fmax increases monotonically with 
Psn/e when i? = 50 fmax increases by a factor of 1.25 
when pBoje is increased to 10^^ cm“^ (and / = 32 GHz) 
from pBo/e = Q (/ = 26 GHz). 

To quantify the possible benefits of lattice alignment. 
Fig. I^a) shows the effect of changing 9 on Ib{Vb). As 9 
increases, the position of the current peak shifts to higher 
Vb- The peak current amplitude, Ipeak, decreases as 9 
increases due to increasing misorientation of the spatial 
parts of the wavefunction, see Fig. 5(b). For example, 
our analysis suggests that the amplitude of the resonant 
peak could be ~ lOx larger for an aligned device. How¬ 
ever, for undoped samples, the PVR increases with in¬ 
creasing 0, see inset in Fig. 5(b), converging to a value 
of 3.4 as 9 approaches 2°: at higher 9, more states are 
available to tunnel resonantly at the current peak m- 
For the doped samples {pbd/& = 10^^ cm ^), however, 
the valley current is close to 0 for all 9, thus the PVR 
is consistently large, see Fig. 5(c). Consequently, the 
increase in current magnitude, which results from align¬ 
ment, leads to higher frequencies without the reduction of 
power that is associated with undoped samples. We find 
that, generally, Rn (oc ^jfmaxi see Eq. (7)) decreases 



FIG. 5. (a) IbfVb) curves calculated for samples with mis¬ 
alignment angles 6 = Q° (black curve), 0.5° (blue curve), 0.9° 
(green curve), 2° (magenta curve), taking pBo/e = 0 cm“^ 
and Nl = 2. (b) Gurrent amplitude at the peak vs misalign¬ 
ment angle, 6. Inset: PVR of IbfVb) vs 9. (c) Ib{Vb) calcu¬ 
lated when 0 = 0° and pBoje. = 10^^ cm“^. Inset: Rn{9) 
for undoped (upper curve) and pbdI& = 10^^ cm“^ (lower 
curve) diodes, (d) fmax vs R curves calculated for an aligned 
sample (black curve) and misaligned sample with 6 = 0.9° 
(green curve), when pBo/f^ = 10^^ cm“^. Gurves are shown 
solid over the range of R presently obtainable in GRTDs and 
dashed for R values that could be achieved by future device 
designs. For all curves, pToje. = 0 cm“^. 


with decreasing 0, Fig. 5(c) inset, and with increasing 
Pbd, meaning that the oscillation frequency is highest 
for 0 = 0° and when pBu = 0. Fig. 5(d) shows that 
perfect alignment could increase fmax by a factor of ^ 2, 
i.e. for i? = 50 H, fmax = 65 GHz when 0 = 0° compared 
to 32 GHz when 0 = 0.9°. The numerical results diverge 
from the small signal analysis power law of fmax oc 
as becomes small, see black curve of Fig. 5(d), and 
it becomes necessary to vary V to induce current oscilla¬ 
tions. 

In conclusion, we have investigated the performance of 
GRTDs as the active element in RLG-based oscillators. 
Our simulations predict that these devices could oscillate 
at mid-GHz frequencies, by careful design of the RLG cir¬ 
cuit. We have also quantified the effect of changing the 
parameters of the GRTD. Reducing the barrier width (a 
modest change to the structure of existing devices) in¬ 
creases the tunnel current, and thus raises the oscillation 
frequency by an order of magnitude. Adjustment of the 
doping of the electrodes can also be used to tune and 
enhance the oscillation frequency. Finally, we have con¬ 
sidered the effect of misalignment of the graphene elec¬ 
trodes: in devices with aligned lattices, frequencies ap¬ 
proaching 1 THz may be attainable. GaAsInAs/AlAs 
RTDs [T] with two layer-thick barriers have similar peak 
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currents and voltages as the GRTD reported here. We 
therefore expect that the GRTD will produce similar EM 
emission power (~10 pW). Our results illustrate the po¬ 
tential of graphene tunnel structures for making the HE 
components in graphene electronics. 
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